Extremal dynamics on complex networks: Analytic solutions 
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The Bak-Sneppen model displaying punctuated equilibria in biological evolution is studied on 
random complex networks. By using the rate equation and the random walk approaches, we obtain 
the analytic solution of the fitness threshold x c to be / + 1), where (k) f = (k 2 ) / (k) (= (k)) in 

the quenched (annealed) updating case, where (k n ) is the n-th moment of the degree distribution. 
Thus, the threshold is zero (finite) for the degree exponent 7 < 3 (7 > 3) for the quenched case 
in the thermodynamic limit. The theoretical value x c fits well to the numerical simulation data in 
the annealed case only. Avalanche size, defined as the duration of successive mutations below the 
threshold, exhibits a critical behavior as its distribution follows a power law, P a (s) ~ s~ 3//2 . 



I. INTRODUCTION 



Punctuated equilibrium is an evolution taking place 
through intermittent bursts of activity separating rel- 
atively long periods of quiescence, which can be often 
found in ecological systems 0, Q|. Bak and Sneppen 
(BS) [3j introduced a simple model to mimic such an 
evolution. The basis of the BS model is to focus on a 
minimal set of variables that capture the basic features 
of punctuated equilibrium while ignoring all other details. 
In the original model, N species are arranged on a one- 
dimensional chain with periodic boundary conditions. A 
fitness value Xi is assigned to each site i (species) on the 
chain, which is a random variable selected in the inter- 
val [0,1]. Evolution in the ecological systems is modeled 
as follows: At each time step, the ecological system is 
updated by locating the site with the lowest fitness and 
mutating it by assigning new random numbers to that 
site and the K — 1 nearest neighboring sites. Subsequent 
updating of the lowest fitness value generates spatial and 
temporal correlations and displays punctuated equilibria. 
A distinct feature arising through this dynamics can be 
found in the distribution of fitness values. After a tran- 
sient period, the distribution of the fitness values has a 
discontinuity at a threshold x c ~ 0.67; its elements are 
zero up to x c and almost the same constant beyond x c . 
The threshold x c is self-organized. 

A mean field version of the BS model was intro- 
duced Q, in which updated are the minimum fitness 
value as well as the fitness values of other K — 1 sites 
selected at random in the system. Such a modified 
model enables one to solve the problem analytically. The 
threshold was obtained to be i c = 1/K in the limit 
N — > 00. Also the notion of avalanche was introduced to 
quantify the correlation between bursts of evolutionary 
activity. Avalanche size is the time interval between two 
successive occasions where no fitness value is less than a 
given value. It was proposed based on the branching pro- 
cess analysis Hand later derived by using the random 
walk approach |(J Q that the avalanche size distribution 
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when the given 



follows a power law as P a ( s ) ~ s 
value is chosen as the threshold x c . 

Ecological systems in real world are complex. Inter- 
actions between individual species are not as simple as 
one-dimensional, but form a complex network. Thus, it 
would be interesting to extend previous studies of the 
BS model performed in the Euclidean space to complex 
networks such as scale free (SF) networks, while it is 
still controversial whether the ecological systems such as 
food webs are SF- networked systems @. SF networks 
mean that the number of connections to each species, 
called degree k i n g raph theory, follows a power law, 
Pd(k) ~ A; -7 0, lldi [Til Il2| . While such an extension 
is natural, only few studies in that direction have been 
performed so far. Christensen et al. Il3l have studied 
the BS model on random networks (14) . Kulkarni et 
al. |l5j studied it on the small- world network introduced 
by Watts and Strogatz ,16|. Moreno and Vazquez 
studied the BS model on SF networks with 7 = 3, ob- 
taining that the threshold x c is given as x c = (k) /(k 2 ) by 
using the heuristic argument similar to the one used in 
the contact process, where (• • • ) denotes the average over 
the degree distribution. They found that x c depends on 
the system size N as x c ~ 1/lnAT, so that it vanishes 
in the limit N — > 00. The iV-dependent behavior was 
obtained numerically at 7 = 3, so that the result may be 
rooted from logarithmic correction. Recently, Lee and 
Kim 0| also studied the same problem on SF networks 
but with general 7 > 2. They obtained that the thresh- 
old is given asx c = ((fc) + l)/((fc + l) 2 ) by using heuristic 
arguments. Thus the threshold vanishes for 7 < 3 and 
finite for 7 > 3. The interesting feature they obtained is 
the crossover behavior in the avalanche size distribution 
between two different power-law behaviors. 

Here we study the BS model on random SF networks 
analytically by using both the rate equation and the ran- 
dom walk approaches [l9| . By random SF networks, we 
mean the SF network with no degree-degree correlation. 
The rate equation is set up for the case that updating 
of the fitness values is carried out not only at the ver- 
tex with the smallest fitness value but also at its nearest 
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neighbors, which is called quenched case. We also com- 
pare the quenched case with the annealed case, where 
updating is carried out at the vertex with the minimum 
fitness values as well as the vertices randomly chosen over 
the entire system and its number is equal to the degree 
of the vertex with the minimum fitness value. It is note- 
worthy that the number of vertices updated is not con- 
stant in complex networks, but depends on the degree 
of the vertex with the minimum fitness value. Thus, the 
analytic approach of the BS model is not as simple as 
the case in the Euclidean space. Here by applying the 
rate equation approach as well as the random walk ap- 
proach, we obtain the fitness threshold analytically to 
be x c = l/((k)f + 1), where (k) f = (k 2 )/(k) in the 
quenched case and (k) / = (k) in the annealed case. The 
avalanche size distribution is obtained to be P a (s) ~ s~ T 
with t = 3/2 at x c for 7 > 3. 

II. RATE EQUATION APPROACH 

In SF networks, it would be essential to take it into 
account the fact that vertices with different degrees ex- 



perience different updating frequencies. Let us denote 
by fk the probability that the vertex with the smallest 
fitness value has degree k. pk{x) is the distribution func- 
tion of fitness values at the vertices with degree k. Here 
we set up the master equation for the fitness distribution 
for the quenched updating, following p| ■ First we define 
a quantity Qk(x), the accumulative distribution of the 
fitness values at vertices with degree k: 

Q k {x) = f dx'p k {x'). (1) 

Thus the fitness distribution pk(x) — — -^Qk{x). Then, 
we have 

fk = - f dx 'jri {Qk(x')} Npk I] {Qk'(x')} NPk ' , (2) 

where pk denotes the degree distribution Pd{k) for sim- 
plicity. Note that ^ fe fk = 1 is satisfied. The evolution 
equation for the fitness distribution at vertices with de- 
gree k is written as 



p k {x,t + 1) = pk(x,t) 



fk 

Np k 
kpk 



£{Q k (x)} Npk i\k^k{Qk'(x)r Pk ' 
fk 



Pk (x, t)-J^ {Qk(x)} Npk Ylk^k {Qk'(x)} 1 
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(k) L 

kpk fk"k" 
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Np k *-f (k) N Pk 
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where the second term of the right-hand side (RHS) of Eq.@ represents the update of the minimum fitness when 
it locates at the vertex with degree k. The third term does the update of the fitness value of a vertex with degree 
k induced by a nearest neighboring vertex with degree k" which has the minimum fitness value in the system. The 
factor kpk/(k) comes from the conditional probability P(k\k") that the vertex with degree k is connected to the one 
with k", which is relevant to the quenched case. In the annealed case, the factor is replaced with pk simply. The last 
two terms do the addition of new fitness values Q. 

The stationary solution in the limit t — > 00 can be solved by using pk(x) — —-^Qk( x ) an d taking the integral over 
[x, 1] of the whole formula as the integral equation, 



l dx '^ {Qkix ' )}NPk n {Q»>(x')} Npk ' 

(k) f kp k 



1 



k(k) f 



[ Np k (k)(N Pk -fk)N 



(k) (N Pk - fk) 



fk + 



(k) 



Np k 
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(4) 



where As done in the threshold x c is determined by the 

comparison of the first term with the second term of Eq. 
^2, ^2, QJ in their absolute magnitudes. To proceed, let us first 

W = kpk anci ( k 'f^£^ k $ k ' ( 5 ) assume that the second term is dominant compared with 

fc=l k=l 
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the first term. Then, we obtain that within the leading 
order, 



Qk(x) 
which leads to 




+ 1 (1-z), 



Pk(x) = 



(6) 



(7) 



This result holds when Q k (x) is less than 1 by more than 
0(l/N). In other words, 



x > x c , k + 0(1/N) 



(k)fk 



(k) f k + (k) f k Pk 



0(1/N), (8) 



where x c ,k is the threshold for a given k. Eq. JJJ) indicates 
that pk(x) does not vanish as N — > oo. We show that 
x C) fc does not depend on k in Appendix. Thus we denote 
x c .k = x c simply. 

For x < x c , the second term of Eq. can be ignored 
and we have 



1 dx ,Np kPk {x') 



Qk(x') 



k(k) 



(k}N 



. | A + 5*^ ,(!-.,. (9) 



Next we use the fact that the integral J^ 1 • • • is very 
small and the interval of the integration is replaced by 
J " . By setting x = 0, we have 



npkpk oc A-. 



(10) 



On the other hand, by differentiating Eq. (JHJ at x = 0, 
we obtain 



n-Pfc/Ofe = ./fc + 



kp k (k) 



f 



(k) 



Combining Eas-HlUfl and l|ll|) . we obtain that 

kp k 



fk 



(11) 



(12) 



which is supported by numerical simulations shown in 
Fig. 1. Note that the result of Ea. H12|) is for the quenched 
case. In the annealed case, similar calculations lead to 
fk = Pk- Thus we obtain (k) , = (k) , = (fc 2 ) / (k) in the 
quenched case and (k)f = (k) ^ a — (k) in the annealed 

case. (k)f t q diverges for 7 < 3 as ~ m 
limit N — > 00. 
We also have 



fc 



pk 



1 



n (k) 



(13) 




FIG. 1: The data of f k (o) and kp k /(k) (x) for the Erdos- 
Renyi network (a) and for the scale-free network with 7 = 3.6 
(b) in the quenched cases. To generate the scale-free network, 
we use the static model ppjl . Both networks have the average 
degree (fc) = 4 and the system size iV = 10 6 . 



Thus, p k - ^(jV-Cr-a)/^-!)), and converges to as 
N — > 00. The convergence rate is slower than the rate 
of 1/N that appears in Euclidean space. Finally, from 
Eq. J5J, the threshold can be expressed simply as 



(k), 



V 



(14) 



which does not depend on k for both updating rules. 
This result is different from the previous results 0j • 
The threshold formula is reproduced by the random walk 
approach in the next section. 



III. RANDOM WALK APPROACH 

The random walk approach was first introduced in 
0, , and it is useful for calculating the avalanche size 
distribution. The threshold can be also obtained. Let 
q\ (t) be the probability of having an avalanche with size 
t, which is defined as the duration of time throughout 
which the minimum fitness value is smaller than a given 
threshold value A. A can be chosen arbitrary. Later we 
find that q\(t) follows a power law when A is equal to the 
threshold x c . The corresponding generating function is 
defined as x( z ) = J2t>o^ii) zt ■ Then x( z ) satisfies the 
self-consistent equation [arill^. 

X (z) = zg( X (z)). (15) 

In the previous study 0, the generating function g(z) 
was given as - A) K "V = (l-\ + Xz) K 

when K fitness values are updated randomly. However, 
in the case of SF networks, the number of vertices up- 
dated at each time step is not constant, but it depends 
on the degree of the vertex with the minimum fitness 
value. In this case, the generating function g[z) is given 
as 



g(z)^J2^-^ 1 ~ x + Xz ^ 



(16) 



fc=i 
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where fk was defined as the probability that the mini- 
mum fitness locates at the vertex with degree k. 

What we do next is to solve q\{t) by using Eos. (1151161) 
To proceed, we use the Lagrange's inversion formula [23 . 
0, 



where z = u>/g(u>), provided that u>/g(u>) is analytic near 
lu = and h(u>) is an infinitely diffcrcntiable function. 
Here we choose h(uo) — to and lo(z) = x( z )- Then, 
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(17) 
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Since x( z ) = q\(t)z l , ?a(£) can be obtained by using the Stirling's formula as 



-. OO OO / t \ 
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(EL, 



-A* -1 (1 — A)^* =1 ki ~ t+1 . 



(18) 



(19) 



Note that fci is the degree of the vertex with the minimum fitness value at updating time i plus 1, which occurs with 
the probability _ i . Thus in the limit t —> oo 



(20) 



Substituting Eq. into Eq. JTHJ) yields 

1 - A 



A 



A(l-A)W/((fc) / + l)W/ +1 
((fc) f ) <fe> ^ 



t~ 3/2 + o{r 5/2 ). 



(21) 



The quantity in the square bracket is 1 when Therefore, the avalanche size distribution behaves as 

A =«7^ (22) 

which is equal to the threshold x c previously obtained 
via the rate equation approach. Then, q Xc (t) ~ t~ 3 / 2 . 
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FIG. 2: (a) and (b). The avalanche size distribution for the 
regular network with the degree distribution of the <5-function. 
In the quenched case (a), the avalanche size distribution fol- 
lows a power law when A is chosen as 0.253 (o), larger than 
the theoretical value l/({fc) + 1) = 0.2 (□). In the annealed 
case (b), the theoretical value x c = 0.2 (o) works well to gen- 
erate the power- law behavior of P a (s). (c) and (d). Same plot 
for the Erdos-Renyi network [r4| . In the quenched case (c), 
the power-law behavior of Pa(s) occurs at x c — 0.207, larger 
than the theoretical value, x c = 1/6. In the annealed case 
(d), the theoretical value :r c =0.2 generates the power-law be- 
havior, (e) and (f). Same plot for the scale-free network with 
7 = 3.6. In the quenched case (e), the power-law behavior of 
P a {s) occurs at x c = 0.1, which is larger than the theoretical 
value x c = 0.08. In the annealed case, the theoretical value 
x c = 0.2 yields the power-law behavior of P a {s). The mean 
degree is fixed to be (k) — 4 and the system size is N — 10 6 
in all cases. The straight lines have slope -3/2 in all cases. 



IV. NUMERICAL RESULTS 

We check the analytical solution of q\(t) numerically 
for several networks. The theoretical formula of x c is 
tested through the criticality of the avalanche size dis- 
tribution. First, the random network with the degree 
distribution Pd(k) — i5fc,fc 7 called the regular network, is 
constructed and the dynamics of the BS model is per- 
formed on that network, ko = 4 is taken for numerical 
simulations. In this reduces to x c = l/(feo + 1) 



simply in both the quenched and annealed cases. In the 
quenched case (Fig.[2Ja)), the avalanche size distribution 
does not follow a power law when A = x c = 0.2. In- 
stead, the power-law behavior appears at a larger value, 
A « 0.253. In the annealed case (Fig. Hfb)), however, 
it follows a power law at A = x c , consistent with the 
theoretical value. 

Second, for Erdos-Renyi (ER) random graph, where 
the degree distribution is a Poisson distribution, the 
theoretical formula reduces to x c « !/((&) + 2) in the 
quenched case, because (k 2 ) = (k) 2 + (k) in the limit 
N — > oo. Numerical simulations are performed in both 
the quenched and the annealed cases. In the quenched 
case (Fig.|2c)), the avalanche size distribution does not 
follow a power law when A is taken as the theoretical 
value, but it does when A « 0.207. In the annealed 
case (Fig. |5Jd)), the avalanche size distribution follows 
a power law at A = l/((fc) + 1), consistent with the the- 
oretical value. 

Next, for SF networks with 7 = 3.6, which is con- 
structed through the static model [2p|. the theoretical 
value x c = l/((fc)/ + l) w 0.08 in the quenched case. 
Note that (k) f » 11.5 is different from (k 2 )/(k) « 7.06 
numerically due to the strong fluctuations arising in the 
large k region (Fig. 1). Again the avalanche size distribu- 
tion does not follow a power law at the theoretical value, 
but does at A w 0.1. In the annealed case, the avalanche 
size distribution follows a power law at x c = l/((k) f + 1). 

The numerical results for the above three networks in- 
dicate that the mean-field theoretical prediction is not 
good for the quenched case, however, is good for the an- 
nealed case instead. This result is attributed to the effect 
of the temporal and spatial correlation between the ver- 
tices with the minimum fitness value at successive time 
steps, which often occur at the nearest neighbors or at 
the same vertex. Such effect was not counted properly in 
the quenched case, and can be neglected in the annealed 
case. 



V. CONCLUSIONS 

We have studied the Bak-Sneppen model on complex 
networks by using the master equation as well as the ran- 
dom walk approaches. The threshold x c is obtained to be 
x c = l/((fc) / + 1), where (•••)/ 1S tne average over the 
minimum fitness vertices. The avalanche size distribu- 
tion follows a power law with the exponent r = 3/2 at the 
critical point. The theoretical prediction of x c was tested 
numerically for the regular network, the ER random net- 
work, and the SF network. For all the networks, the 
theoretical predictions of x c are in disagreement (agree- 
ment) with the numerical results for the quenched (an- 
nealed) case. The discrepancy in the quenched case is 
attributed to the effect of the temporal and spatial cor- 
relation between the vertices with the minimum fitness 
values at successive time steps. Nevertheless, the for- 
mula x c = (k)/{(k 2 ) + (k)) is newly derived here for the 
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quenched case. Thus when 2 < 7 < 3, x c — > in the 
thermodynamic limit in the quenched case. 



APPENDIX A 

Here we show that x c k does not depend on k. To 
proceed, we suppose Xk ltC < Xk 2 ,c for a certain pair of k\ 
and hi. Then, for a given xq in the range Xk ltC < %o < 
x k 2 ,c, we have with Eq. 



,Np kl p kl (x') 



dx 



< — 7^ Qfci(zo), (Al) 



and 



dec 



,Np k2 p k2 (x') 



<fc) 



1 - 



fe 2 (fc> / 

(fc) TV 



(&) f k 2 Pk 2 
(k) 



Q k2 (x ). (A2) 



Since Q kl (x) and Qfc 2 (x) are of the same order 
based on Eq. ©, the RHS of Eq. l|Alj) and that of 
Eq. (jSU are of the same order. Then we must have 
Pki(x') <C pk 2 (x') for xq < x' < 1, which contradicts 
Eq. J7J) for a; > Xfe 2 . c (> a;o,Xfc liC ). Thus we set x c = Xk, c 
for all k. 
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